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The present work reports the analysis of a possible relationship due to stress transfer between the two 
earthquakes that hit the province of Van, Eastern Turkey, on October 23, 2011 (Mw — 7.2) and on 
November 9, 201 1 (Mw =5.6). The surface displacement field of the mainshockhas been obtained through a 
combined data set made up of differential interferograms from COSMO-SkyMed and ENVISAT satellites, 
integrated with continuous GPS recordings from the Turkish TUSAGA-AKTIF network. This allowed us to 
retrieve the geometry and the slip distribution of the seismic source and to compute the Coulomb Failure 
Function (CFF) variation on the aftershock plane, in order to assess a possible causal relationship between 
the two events. Our results show that the November 9 earthquake could have been triggered by the October 
23 shock, with transferred stress values largely exceeding 1 bar. 

The occurrence of a seismic event perturbs the stress field within the Earth crust, altering the probability of 
occurrence of a second strong event along nearby faults. This phenomenon is known as fault interaction, and 
the physical mechanism at its basis is known as the Coulomb Stress Triggering'. In literature it has been 
shown that even weak increments of stress can trigger large seismic events^, if the relative seismogenic structures 
are close to the end of their seismic cycle; furthermore, the interaction between two seismic events may occur both 
within long interseismic periods'* and smaller time intervals'* ''. Stress field perturbations cannot be instrumentally 
detected, but can be numerically simulated by means of mathematical models^'^ Several studies demonstrate 
that the distribution of aftershocks and the variations of the seismicity rates are generally explained in terms of the 
Coulomb faUure criterion''^ '" which allows to interpret a shear stress exceeding a certain fraction of the normal 
stress as an increased probability of failure'". Several physical-static prediction models based on the Coulomb 
criterion have been proposed" '''; these models have been widely applied for the redefinition of the regional 
seismic hazard following earthquakes'^''"""*. 

The province of Van (Eastern Turkey) has been struck by two strong earthquakes in a few days. The first 
earthquake occurred on October 23, 2011 (Mw = 7.2), 30 km N of Van city, while the second took place on 
November 9, 2011 (Mw = 5.6) few kilometers to the South'"* (Fig. la), in the Edremit subprovince. A significant 
aftershock activity was recorded after the first shock over the Van region (Fig. la). The Mw 7.2 Van earthquake 
occurred along a ENE-WSW fault, previously recognized by Kein^", who attributes to it a main right lateral strike 
sUp mechanism. However the CMT (Centroid Moment Tensor) solution from USGS and GFZ indicates a pure 
reverse fault mechanism (Fig. la)^''^'', while the Harvard CMT solution is a mixed reverse-right lateral mech- 
anism. Concerning the November 9 (hereafter Edremit- Van) earthquake, the Kandilli Observatory and 
Earthquake Research Institute (KOERI) CMT solution indicates a dominant strike slip mechanism^'. 

Turkey is surrounded by three main plates, the African, Eurasian and Arabian ones, characterized by relevant 
tectonic activity. Other two minor plates are also present in this region, the Aegean and the Anatolian ones. The 
relative motion amongst the above mentioned plates has generated some of the major tectonic features of Turkey: 
the Aegean Arc, the West Anatolian Graben Complexes, the North Anatolian Fault Zone, the East Anatolian 
Fault Zone, the North East Anatolian Fault Zone, the Bitlis Thrust Zone, and the Caucasus (Fig. Ib)^*. The study 
area is located where the Arabian Plate, that moves towards north-northeast, collides with the Eurasia along the 
Bitlis Thrust Zone, a complex coUisional boundary located north of the fold-and-thrust belt of the Arabian 
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platform. This suture zone corresponds with the Van Lake region 
from which the other smaller plates move symmetrically away, 
forced by the Arabian plate movement (see Fig. Ib)^^. In this sector 
the GPS velocities indicate a NW-oriented motion of 18 ± 2 mm/yr 
relative to the Eurasia^'", causing an intense seismic activity. This 
latter is associated with ENE-WSW and NW-SE conjugate strike- 
slip faults, of dextral and sinistral character, parallel to the North and 
the East Anatolian fault zones. Among the major structures in the 
area the (^aldiran Fault, Erci^ Fault, Igdir Fault, Malazgirt Fault, 
Siiphan Fault, Kagizman Fault Zone, Tutak Fault Zone and 
Northeast Anatolian Fault Zone can also be mentioned. These struc- 
tures are compatible with NNE-SSW compression and WNW-ESE 
extension in the region that well-matched the reverse mechanism 
related to the Van earthquake. 

Several large earthquakes have struck this region both in historical 
and instrumental ages. The most relevant known that damaged the 
cityofVan was the 1648 event (March 31, Mw ~ 6.6), which has been 
associated with the E-W trending Gurpinar thrust, located south of 
Van^^. The 1715 event (March 8, Mw = 6.6) was located between the 
eastern termination of Derik fault and Van^", whilst the 1903 (April 
28, Mw = 7) Malazgirt earthquake occurred on a NNE-SSW trend- 
ing structure'''. One of the most recent and destructive event was the 
November 24, 1976 (Ms = 7.3), Caldiran earthquake, causing more 
than 4,000 casualties and intense damaging in a 2,000 km^ area'". 
The responsible fault was a N 110°-135°, 50 km long, right lateral 
strike-slip mechanism, with observed horizontal displacements vary- 
ing between 2.5 and 3.5 m and a vertical displacement of about 
50 cm. 

We have applied Differential Interferometric SAR (DInSAR) tech- 
nique to investigate the surface displacement due to the 20 11 Mw 7.2 
and Mw 5.6 shocks. Different datasets have been used. In order to 
study the mainshock, a pair of COSMO-SkyMed Stripmap data has 
been processed, together with two ENVISAT SAR images. An addi- 
tional dataset composed of two TerraSAR-X pairs was available to 
measure the surface effects induced by the second event and to infer 
its fault geometry. Finally, in order to remove the topographic phase 
contribution, we used the Shuttle Radar Topographic Mission 
(SRTM) digital elevation model. To better constrain the deformation 
field, we integrated in our analysis the DInSAR dataset with coseis- 
mic displacement field from 20 GPS stations belonging to the 
Turkish CGPS network, whose static offsets was assembled by Dr. 
Rahsan Cakmak (TUBITAK MRC EMSI) and are available at http:// 
supersites.earthobservations.org/van.php. 

Results 

DinSAR processing has been applied to a SAR dataset composed of 
three pairs of Very High Resolution (VHR) X-band images (in strip- 
map mode" "), one acquired from COSMO-SkyMed (10/10/2011- 
23/10/2011), i.e. related to the main event, and the other two from 
TerraSAR-X (9/11/2011-20/11/2011 and 31/10/2011-11/11/2011) 
satellites. The frame available from COSMO-SkyMed partially cov- 
ers the epicentral region, with the epicenter of October 23 mainshock 
located few km out of the area of measurement (Fig. 2). The mea- 
sured surface movements have a maximum deformation along the 
LOS (Line Of Sight) of about 1.0 m in the uplifting sector, while in 
the southern portion of the interferogram a subsidence up to 0.02 m 
has been detected. As the COSMO-SkyMed images did not cover the 
northern section of the deformed area, a pair of SAR images acquired 
by ENVISAT platform (22/07/2011-19/11/2011) has been used in 
order to better assess the induced surface displacement in the upper 
region of the displacement field, (Fig. 2). Despite the ENVISAT data 
pair covers a time interval including the November 9 event, the 
northern portion of the displacement field can be considered free 
of the effects of the second earthquake because of to the large epi- 
central distance. Two TerraSAR-X interferograms, along descending 
and ascending orbits (Supplementary material Fig. SI), have been 



used to investigate the effects of the November 9 earthquake. These 
interferograms clearly highlighted that the November 9 event did not 
contribute to the deformation induced by the main shock on October 
23. 

TerraSAR-X interferograms show different fringe patterns, prob- 
ably due to the strike slip mechanism, that implies most of the move- 
ment along an horizontal axis. Unfortunately TerraSAR-X data do 
not allow to solve the ambiguity of the E-W and N-S conjugate 
planes, either because the coseismic displacement field is partially 
detected, or due to the moderate magnitude of the buried fault. 
Nevertheless, the E-W plane solution can be considered the most 
probable as suggested by Akinci et al.'''' and supported by field 
evidence'"'. 

We retrieved the source model for the Van earthquake from a two- 
step joint inversion of DInSAR and GPS data. In a first step, the 
geometry and extension of the source have been determined assum- 
ing uniform slip on the fault plane. The probability distributions of 
the fault plane parameters (Fig. S2) have been obtained by evaluating 
the Bayesian integrals with a numerical Markov-Chain Monte Carlo 
integration scheme'^. The most probable model (Table 1) has a strike 
of 252 degrees and dips NNW at 50 degrees. The strike value is 
consistent with the orientation of the CMT focal plane and with 
geodetic models'"" '^. The modeled dip angle is slightly larger than 
the value from Harvard CMT, but is consistent with geodetic indica- 
tions of a steeply dipping fault in the range 40-55 degrees'"'. The 
geodetic magnitude for the uniform-slip model is M = 7.08 if a shear 
crustal modulus of 30 GPa is assumed. 

After the definition of the fault geometry, a linear inversion has 
been performed in order to estimate the slip distribution on the 
rupture plane. In this step, the size of the fault plane has been 
extended in order to better account for slip heterogeneities. The 
resulting slip distribution is shown in Figure 3 and has a cumulative 
geodetic magnitude of M = 7.12, in good agreement with the seis- 
mological estimates. The resulting pattern shows a deep high-slip 
patch, with peak slip (3.8 m) located at 15 km depth, with the epi- 
central location (provided by KOERI) located at its lower boundary; 
the location and extents of the bulk slip area is consistent with results 
by Fielding et al."". This geometry suggests that the rupture nucleated 
in the deep central portion of the fault and migrated towards the 
surface. The direction of slip vectors show dominant reverse-left 
mechanism; a patch with moderate reverse-right slip values at the 
southeastern tip of the modeled fault is probably an artifact, because 
this edge of the fault is beneath water and therefore is not constrained 
by SAR data. The modeled deformation field (Figures 4a and 4b) 
reproduces correctly the observed displacement; residuals of a few 
centimeters found on the eastern portion of the fault plane surface 
trace may be ascribed to a change in the strike direction''^, and/or to 
the activation of secondary shallow ruptures'"*. 

Once defined the source model for the Van earthquake, the role of 
this event in promoting the rupture of the Van-Edremit earthquake 
through a Coulomb Failure Function (CFF) analysis can assessed. 
The CFF is evaluated by computing the incremental stress tensor 
produced by the elastic dislocation of the Van earthquake, projecting 
it on the rupture plane of the Van-Edremit earthquake and evalu- 
ating the relative contributions of the normal and shear stresses. 
Positive or negative variations of the CFF indicate that the perturba- 
tion to stress field is acting to promote or oppose the rupture, respect- 
ively. Since there is no information available to reliably identify the 
rupture geometry of the Edremit-Van event between the two 
Harvard CMT conjugate planes, we computed the CFF variation 
on both planes. The two planes are centered on the hypocentral 
location given by KOERI (43.234N, 38.430E), with dimensions 15 
X 10 km^, covering a depth range of about 0-10 km. In Fig. 5 we 
show the variation of CFF resulting on the two CMT planes. Both 
planes are loaded with positive CFF variations, with average values of 
1.3 bar and 1.1 bar for the E-W and N-S planes, respectively. Peak 



SCIENTIFIC REPORTS | 4:3959 | DDI: 1 0. 1 038/srep03959 



2 




43°0'0"E 



Eurasic 



asian Plate 



Black Sea 



44°0'0"E 




4 . 

Arabian Plate W 




200 
I 



Figure 1 | (A) Aftershocks distribution following the October 23, 2011, mainshock: in blue up to November 9, 2011, earthquake; in yellow since 
November 9 and till November 30, 20 1 1 . Red lines represents the main active faults in the area, modified from SHARE "Seismic Hazard Harmonization in 
Europe" (www.share-eu.org). Figure lA is not released under a Creative Commons Attribution-NonCommercial-ShareALike 3.0 Unported License. This 
image is licensed under a separate, Creative Commons Attribution-ShareAlike 3.0 Unported License. To view a copy of this licence visit http:// 
creativecommons.org/licenses/by-sa/3.0/ (B) Tectonic setting of the Eastern Mediterranean region. White and black arrows indicate the plate motion 
versus. Maps created with ARCGIS 10 software. 
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Figure 2 | Mosaic of unwrapped differential interferograms relative to the October 23, 2011 earthquake, from COSMO-SkyMed and ENVISAT 
coseismic pairs (south and north of the dashed line respectively). The ENVISAT dataset has been scaled to compensate the small discontinuity with 
respect to COSMO-SkyMed data mainly due to a different reference point selected for the unwrapping procedure, and the different LOS of the two 
satellites. The displacement is in LOS geometry. Maps created with ARCGIS 10 software. 



CFF values (2.6 bar for the E-W plane and 3.0 bar for the N-S plane) 
are found on the deep portion for both planes. 

Discussion 

We have investigated the Van earthquake (October 23, 2011) by 
applying DInSAR technique and exploiting the capabilities of the 



Table 1 [ Fault parameters for the October 23, 201 1 earthquake 
resulting from the joint inversion of the DInSAR (COSMO-SkyMed 
and ENVISAT) and GPS dataset 

Most probable model 95% confidence interval 


UTM X (km) 


357.7 


-4.0 




+4.0 


UTM Y (km) 


4640.9 


-4.0 




-H3.0 


Top depth (km) 


9.8 


-2.7 






4-3.3 


Strike 


252.0 


-7.0 






-H6.0 


Dip 


50.3 


-5.5 






+4.5 


Rake 


74.4 


-8.8 






-H8.0 


Slip (m) 


3.8 


-1.2 






+ 3.4 


Fault length (km) 


20.7 


-7.0 




+ 7.7 


Fault width (km) 


16.9 


-9.3 




+ 9.8 



VHR SAR data acquired by COSMO-SkyMed and the C-band 
ENVISAT images. The COSMO-SkyMed differential interferogram 
has allowed to retrieve the surface displacement field. Since the 
DInSAR data does not cover the region located to the north of the 
event, we integrated SAR measurements with coseismic GPS offsets 
from the Turkish geodetic network. 

The inference of the October 23 seismic source on the basis of 
geodetic data has been the subject of a set of studies"""'". Fielding 
et al.^' modeled the rupture as a single fault plane using SAR, GPS 
and seismic waveform data. On the other hand, Elliott et al.^* mod- 
eled surface displacements from COSMO-SkyMed and ENVISAT 
with a pair en-echelon fault planes. Both studies indicate a deep 
rupture, with significant slip at 8 km depth, a dip range between 
40 and 54 degrees. In our analysis, we modeled the fault with a single 
plane. Our most probable solution has a mixed thrust-left slip on a 
plane dipping NNW at about 50 degrees, consistently with the cited 
results even if slightly larger than the dip from Harvard CMT. 
Moreover the strike is consistent with seismological and geodetic 
estimate. Most of the slip occurs on a patch of about 15 X 10 km^, 
approximately located at depths between 10 and 18 km. The location 
and extents of the high slip patch is consistent with the peak slip area 
obtained by Fielding et al.". Probability density functions for model 
parameters (Fig. S2) shows that the position, geometry and along- 
strike extents of the fault plane are well resolved, while trade-offs 
exist among depth, slip and along-dip extents, which is a well-known 
limitation of geodetic source inversions''. The highest residuals 
(Fig. 4a) are in the central part of the northern edge of the 
COSMO-SkyMed frame, where boundary effects maybe present, 
and along the eastern portion of the fault surface trace, supporting 
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the hypothesis of a discontinuity in the strike direction as suggested 
by Elliot et al.^^. GPS horizontal offsets (Fig. 4b) are well recovered by 
the model, while the largest offsets (at sites MURA and OZAL) in the 
vertical component are overestimated, even though residuals are well 
within a standard deviation. 

Finally, we investigated the effect of the perturbation to regional 
stress field induced by the mainshock on the November 9 earthquake 
through the evaluation of the CFF variation. For the Edremit-Van 
earthquake, Harvard CMT solution provides a N-S and an E-W 
conjugate planes, both with strike-slip geometries. The E-W plane 
is the most convincing; this conclusion is supported by field evi- 
dences found by Selcuk et al.^*. However, since it is not possible at 
present to reliably solve the ambiguity between the two geometries, 
the CFF variation has been computed on the two conjugate planes 
(Fig. 5). We found that both are loaded with stress levels up to 2.5- 
3.0 bar, largely exceeding the threshold value of 0.1 bar that is widely 
used to assess effective triggering of seismic events^'^'*. These findings, 
in agreement with previous results by Akinci & Antonioli^^, support 
the hypothesis that the Edremit-Van earthquake has been actively 
promoted by the October 23 Van event. 

Methods 

In order to study the surface effect of the two main events, we have applied DInSAR 
method to a dataset composed of four pairs of images. The first, one was acquired by 
the Italian COSMO-SkyMed satellite, the second one by the European ENVISAT 
platform, and the two other pairs by the German TerraSAR-X. COSMO-SkyMed and 
TerraSAR-X satellites are equipped with an active X-band microwave sensor, whUe 
ENVISAT SAR is a C-band system, all capable to achieve cloud-free and day-and- 
night land observations. The X-band images are acquired in Stripmap mode at 3 m 
spatial resolution, and the ENVISAT data are in IS6 mode {https://earth.esa.int/ 
handbooks/asar/), at about 20 m spatial resolution. The COSMO-SkyMed and 
ENVISAT data, both along ascending orbit, have been used to investigate the October 
23 mainshock. In particular, the post-event image of COSMO-SkyMed is dated 
October 23, only few hours after the seismic event, while pre-seismic one was acquired 



on October 10, 2011. The ENVISAT scenes were taken on July 22, 2011, and 
November 19, 201 1. It is worth to note that the time interval between these two images 
includes also the strong aftershock occurred on November 9. We selected from 
ENVISAT interferogram the portion that complete to the north the surface dis- 
placement measured from COSMO-SkyMed dataset (see Fig. 2), and which was not 
affected by the November aftershock, as confirmed by TerraSAR-X interferograms. 

TerraSAR-X images have been used to measure the deformation caused by the 
November 9, 2011, earthquake. In this case for the descending pair, the pre-seismic 
image is acquired few hours before the earthquake, while the post-seismic one is dated 
November 20, 2011. Concerning the ascending pair, the pre-seismic image is ten days 
before the earthquake (31/10/2011), the post-seismic image two days after (11/11/ 
2011). 

In order to improve the signal-to-noise ratio, a 2 by 2 multi-look factor in slant- 
range and azimuth has been applied to the COSMO-SkyMed and TerraSAR-X 
interferograms, with a square pixel of about 6 m. The 90 m Shuttle Radar 
Topography Mission (SRTM) digital elevation model has been used to remove the 
topographic contribution of the interferometric phase. Since the temporal baseline is 
sufficiently short, both interferograms maintain a good coherence, allowing capturing 
most of the coseismic deformation. Before the retrieval of the Line Of Sight (LOS) 
displacements using a minimum cost flow phase unwrapping algorithm, the phase 
noise has been reduced applying an adaptive filter^. 

Concerning the lower resolution data coming from ENVISAT, we computed the 
interferogram by using a multi-look factor equal to 1 X 5, in range and azimuth 
respectively, and by applying the same phase noise reduction and unwrapping 
methods of X-band data. DORIS precise orbital data have been used to correct 
possible orbital fringes. Despite the low interferometric coherence of ENVISAT 
interferogram, caused by the long time interval (temporal baseline), the unwrapped 
data have contributed to improve the spatial coverage and to better constrain the 
model retrieval. 

We jointly inverted DInSAR (COSMO-SkyMed and ENVISAT) and GPS data 
using a two-step method: we first inferred the geometric features of the Van earth- 
quake source with a non-linear inversion assuming a uniform slip, then we retrieved 
the slip distribution on the fault plane with a linear inversion. For the inversion 
procedure, a combined DInSAR deformation field has been obtained by selecting the 
ENVISAT field laying in the region not covered by COSMO-SkyMed data and 
merging this subfield with the full COSMO-SkyMed field. In regions where both 
datasets are available, we used only the COSMO-SkyMed dataset that, because of its 
shorter temporal baseline, it is less affected by post-seismic effects. The resulting 
combined dataset has been downsampled by a factor of 0.5 for computational reasons. 
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Figure 4 | (A) Observed (left) and modeled (center) DInSAR deformation fields. Residuals (right) are defined as the difference between observed and 
modeled displacements. The surface projection of the source model is also shown as a dashed box. The vertical shift between the COSMO-SkyMed and 
ENVISAT displacement is due to the different absolute references for the two datasets. (B) Observed and modeled horizontal (left panel) and vertical 
(right panel) GPS deformation fields. Error ellipses in the left panel and vertical bars in the right panel correspond to 68% confidence levels. The GPS 
offsets show a compressional deformation regime, consistently with the expected thrust mechanism; the largest offset of about 6 cm is found at the near- 
field site MURA. Maps created with the GMT software. 



For the nonlinear inversion procedure we modeled the deformation field for a 
given fault plane using the analytical expressions from Okada'*\ assuming an 
homogeneous Poisson half-space. To take into account the different characteristics of 
DInSAR and GPS datasets, we applied different weights to each in order to obtain 
similar contributions to the total misfit for models that give an equivalent fit to the 
data. We used a Bayesian inference scheme to estimate the probabUty distribution 
functions (PDF) of fault parameters through a Markov-chain Monte Carlo integ- 
ration'=. Since the COSMO-SkyMed and ENVISAT datasets have different absolute 
references, we included in the model parameters a "bias" value between the two 
datasets. 



Once the geometry is fixed, we estimated the slip distribution on the fault plane 
with a linear inversion. We extended the fault plane resulting from the uniform 
inversion up to 40 X 30 km"* in order to better account for (eventual) slip hetero- 
geneities, and subdivided it into patches of size of about 1.5 X 1.5 km"* We composed 
a Green Function matrix by imposing a unitary slip on each patch and computing the 
corresponding deformation fields for DInSAR and GPS data, according to the Okada 
analytical expressions'^^ The slip distribution is then recovered by simultaneously 
minimizing (in the least-squares sense) the total chi-squared and a discrete 
approximation of the Laplacian to avoid large, unphysical oscillations in slip values. 
The solution is computed assuming a positivity constraint on the model, with an 
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Figure 5 | CFF variations induced by the October 23, 201 1 event on the November 9, 201 1 rupture plane for both Harvard CMT solutions. In the left 
panel, a red box marks the position of the modeled fault plane for the mainshock, while black boxes mark the positions of both CMT nodal planes. Maps 
created with the GMT software. 



optimized version of the Lawson-Hanson method^ *^. The inclusion of a smoothing 
term to damp unphysical model oscillations introduces a tradeoff between data fit and 
solution roughness. We determined the weight of the smoothing constraint by 
evaluating a misfit-roughness tradeoff curve (figure S3). 

We estimated the spatial resolution for the slip distribution on fault plane by 
defining a checkerboard synthetic slip model. The resulting slip distributions have 
been computed with the coverage of SAR and GPS data and inverted with the same 
procedure used for the real datasets. Results (Figure S4) show that a resolution of 
about 5 km can be expected up to —5 km depth, while at greater depths model 
resolution decreases to over 10 km. 

We then investigated the role of the Van earthquake in promoting the rupture of 
the Van-Edremit earthquake by evaluating the Coulomb Failure Function (CFF). 
Using the source model obtained for the Van earthquake, the elastic strain tensor 
corresponding to the seismic dislocation is computed with the analytical solutions 
provided by Okada^^ Using standard relations from elasticity theory, the strain field is 
converted into an incremental stress tensor that acts as a perturbation of the pre- 
existing (unknown) regional stress field. The effect of the perturbation to the stress 
field on a given fault mechanism is then assessed by computing the CFF variation, 
defined as ACFF — At + )J.(Aa„ + Ap), where At and Aa^ are respectively the shear 
and normal incremental stresses, |j, is the friction coefficient and Ap is the pore 
pressure change'**. It is convenient to rewrite this relation as ACFF ~ At + iigAcin' 
where ]ie is an effective friction coefficient taking into account static friction, 
hydrostatic pressure and pore fluid pressure^'^. We assumed ~ 0.4, a value con- 
sistent with laboratory evidences on friction and moderate pore pressure in condi- 
tions where fluids are not fuUy expelled*^. Knowing the values of ACFF on a given 
fault mechanism allows to establish whether the stress field is acting to promote 
(ACFF > 0) or oppose (ACFF < 0) the rupture. As a general rule, a CFF increase of 
0.1 bar (corresponding to tidal load) is considered in literature as effective for 
earthquake triggering^-^". 
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